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Computer simulation of the hopping charge transport in disordered organic materials has been 
carried out explicitly taking into account charge-charge interactions. This approach provides a pos- 
sibility to take into account dynamic correlations that are neglected by more traditional approaches 
like mean field theory. It was found that the effect of interaction is no less significant than the 
usually considered effect of filling of deep states by non-interacting carriers. It was found too that 
carrier mobility generally increases with the increase of carrier density, but the effect of interaction is 
opposite for two models of disordered organic materials: for the non-correlated random distribution 
of energies with Gaussian DOS mobility decreases with the increase of the interaction strength, while 
for the model with long range correlated disorder mobility increases with the increase of interaction 
strength. 

PACS numbers: 72.20.Ee, 72.80.Le, 72.80.Ng 



INTRODUCTION 



MODEL 



Hopping transport in disordered organic materials has 
been extensively studied for the case of low density of 
carriers, but our understanding of charge transport for 
the case of high carrier density is not adequate; theoreti- 
cal studies of the effect of carrier density are scarce [l|, . 
High density of carriers could affect drift mobility in op- 
posite ways. Small fraction of carriers could occupy deep 
states, thus providing a possibility for remaining carriers 
to avoid trapping and acquire much higher mobility. At 
the same time, charge-charge interactions could provide 
an additional energetic disorder in the material. This is 
indeed the case for the simplest model where all charges 
except one are immovable [3|. In this case, the greater is 
the density of static charges (i.e. energetic disorder), the 
smaller is the mobility. 

Usually, in theoretical studies the mean field approxi- 
mation has been used and charge-charge interaction has 
been totally neglected P, 0] . This means that these stud- 
ies dealt only with the effect of filling of deep states: it 
is assumed that the effects of interaction could be later 
effectively included via the mobility dependence on the 
mean local electric field ^-Eioc^, which in turn is con- 
nected to the mean local charge density p by the Poisson 
equation 
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where e is a dielectric constant of the material. This 
line of reasoning totally neglects dynamic correlations. 
In addition, quite frequently true quasi-equilibrium mo- 
bility is formed by the averaging over large domains of 
the disordered material (see, for example, Ref. |^). In 
this situation the very conception of a local (but uniform 
in space) mobility is invalid. For this reason in this pa- 
per we describe a direct dynamic simulation of hopping 
monopolar transport of interacting carriers. 



General description 

For the model of interacting carriers the total site en- 
ergy of a carrier is a sum of a static random contribution 
due to the static intrinsic energetic disorder and fluctu- 
ating contribution due to the interaction with all other 
hopping charges. Tremendous difficulty of the simula- 
tion is a necessity to recalculate site energies and hopping 
probabilities after every hop. In addition, if we consider 
moderate or high density of charge carriers, then they 
affect the distribution of average electric field, which in 
turn leads to the non-uniform average charge density. In 
this situation we cannot use in simulations a finite basic 
sample with periodic boundary conditions. The length of 
the basic sample across the device must be equal to the 
thickness of the transport layer. This inevitably means 
that we have to include charge injection in our simulation 
and analysis of the simulation data becomes much more 
complicated. 

For these reasons we undertook simulation for the sim- 
plified model of a transport layer: we assumed that 
a static compensating charge of the opposite sign is 
uniformly distributed in the transport layer, with the 
density being equal to the average density of carriers. In 
this model the mean local electric field and carrier density 

are uniform in space and ( E\oc ) is equal to the applied 



field E. Hence, simulation for the moderate basic sample 
with periodic boundary conditions could be carried out. 
This particular model is a variant of the "jelly model" , 
very popular for a study of effects of charge-charge inter- 
actions in metals. 

In order to minimize the necessary runtime even more 
we performed simulations for rather high tempera- 
ture kT/a = 0.3 (here is a variance of static intrinsic 
disorder; typically, in organics a ~ 0.05 — 0.1 eV 4]) 
and took into account only hops to nearest neighbors of 
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every site (we considered a simple cubic lattice) . Site en- 
ergies and hopping probabilities have been recalculated 
only for nearest neighbors of sites occupied with carriers. 
This assumption seems to be very good approximation 
for disordered organics where aa ~ 5 — 10 [51], here a ~ 1 
nm is a distance between nearest neighbors (lattice scale 
in a lattice model of a material) and a is an inverse lo- 
calization radius of a wave function of a transport site. 
Miller-Abrahams (MA) hopping rate Q has been used 
for all simulations. In the calculation of electrostatic site 
energies a well-established Ewald approach has been used 

Simulations for two models of a disordered organic 
material have been carried out: the Gaussian Disorder 
Model (GDM) where the static intrinsic energetic disor- 
der in the material is represented by a spatially uncor- 
related distribution of random energies with the Gaus- 
sian DOS [H] and the model of dipolar glass (DG) where 
site energies are strongly correlated 3 • The later model 
is better suited to describe organic materials because 
contributions of static electrostatic sources (of dipolar 
or quadrupolar nature) to site energies are inevitably 
strongly spatially correlated due to the long range nature 
of electrostatic interaction in organic materials jol. [loj . 

Majority of simulations have been performed for the 
basic sample of a simple cubic lattice with the size L — 
40a. Several checks have been made to test an accuracy 
of the simulation by using the basic sample with size 
L — 80a. Even for the DG model and rather 
weak field eaE/a w 0.1 (these are the situations where 
mobility is very sensitive to the size of a sample) a good 
agreement has been found with the case L — 40a. 

Simulation for a particular set of relevant parameters 
was started with some arbitrary distribution of carriers 
and carried out until a stationary state with constant 
mean carrier velocity has been reached. We have checked 
that the particular initial distribution of carriers affects 
only details of the relaxation process but not the value 
of the asymptotic velocity. 

Waiting time simulation 

General features of the simulation procedure are very 
similar to those described in Refs. apart from the 

modification of the waiting time calculation. For the only 
carrier a waiting time r before the next hop is calculated 
as 

T = -^log7, (2) 

where 7 is a random number with the uniform distribu- 
tion in [0,1] and T is the total hopping rate for the site 
where carrier is located at the moment [H. For many 
interacting carriers this procedure has to be modified be- 
cause for any particular carrier, waiting for a hop, the 



total hopping rate now depends on time, reflecting hop- 
ping of other carriers. 

A necessary modification could be introduced in the 
following way. At the start, independent 7„ have been 
generated for all carriers, and the set of Tn has been calcu- 
lated using Eq. ((2]), as well as the set of initial residuals 
Rn = — log 7„ . Then the carrier with smallest waiting 
time Tmin made a hop, its new position has been found 
by the usual way Q , and the time counter was advanced 
by Tinin- For this carrier a new random number "f^™ was 
generated. For all other carriers new residuals 

Rn (3) 

have been calculated and stored (for the hopped carrier 
a new residual is just i?""™ — — log 7"™). Then all site 
energies and hopping rates have been recalculated, the 
new set of waiting times have been calculated 

pncw 

now _ fA\ 
n pncw' ^ ' 

and the whole procedure has been repeated. 



DANGERS OF MEAN FIELD APPROXIMATION 

To illustrate possible dangers of the mean field approx- 
imation we provide results of the transport simulation for 
the simplest case of no intrinsic static disorder (carriers, 
hopping in the empty lattice). In this case the mean field 
solution for the carrier occupation fraction at every site is 
just a constant equal to the average occupation fraction 
p. For the MA hopping rate an average carrier velocity 
for the case of simple cubic lattice and hopping to nearest 
neighbors only is 

WMF = Wo [1 - exp (-ea£;//cT)] (1 -p), wo = oFq, (5) 

where two terms in brackets represent contributions of 
forward and backward hops, factor 1 — p takes into ac- 
count occupation of the final site, and Fq is a hopping 
frequency. Remarkable property of Eq. ([5]) is its indepen- 
dence on the interaction strength. Indeed, in our model 
in the absence of static intrinsic disorder an average total 
charge density is exactly zero everywhere, thus all inter- 
actions enter into consideration only because of dynamic 
correlations. Results of the simulation are shown in Fig. 
[1] For strong interaction deviation with Eq. ([5]) is signif- 
icant. Note also, that for e/ea^E = agreement between 
Eq. ([5]) and simulation data is excellent due to the small 
average occupation fraction p <C 1 used in simulation; 
here dynamic correlations are not important. 
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FIG. 1: Temperature dependence of the mean carrier velocity 
V for the MA hopping rate in the case of no intrinsic disorder 
(dynamic simulation, points). Upper solid line is the mean 
field result. Numbers near the curves show effective charge- 
charge interaction parameter e/ea^E. Average occupation 
fraction is p = 0.047. 



FIG. 2: Dependence of the carrier mobility on the effective 
strength of the charge-charge interaction Uc = /eacr for 
L = 40a, p = 0.008, kT/cr = 0.3, and eaE/a = 0.1. Black 
squares show the data for the DG model and empty circles 
show results for the GDM. For a typical disordered organic 
material Uc — 5. Here ^(0) is the mobility for Uc = 0. 



RESULTS OF DYNAMIC SIMULATION 
General effect of interaction 

We found that in our model for not very high occu- 
pation fractions p < 0.1 carrier drift mobility increases 
with p, exactly as in the case of non-interacting carriers 
[H, Yet modification provided by the interaction is 
still significant (see Fig. [5]). There is a striking difference 
between the effects of interaction (i.e., carrier-carrier re- 
pulsion) on the mobihty in the DG model and GDM. In 
the DG model repulsion between carriers makes mobility 
even greater than in the case of no interaction, while for 
the GDM the opposite situation takes place. This differ- 
ence is not surprising. Indeed, in the DG charge trans- 
port is significantly affected by carrier trapping in deep 
and broad valleys of the energy landscape (good quali- 
tative description is provided in Ref. [ll|). If a carrier 
is trapped by a valley, then the whole valley with many 
sites having low energies becomes blocked for other car- 
riers because of repulsion. Thus, filling of the deep states 
is much more effective in correlated landscape if carrier 
repulsion is taken into account. This is the reason for 
the increase of carrier mobility in DG. No such effect 
takes place for the GDM, and here, evidently, the ef- 
fect of charge-induced energetic disorder is responsible for 
the decrease of mobility in comparison to non-interacting 
case. 



This particular result disagrees with the result of a re- 
cent paper by Zhou et al p^]. They found that in the 
GDM carrier interaction enhances mobility in compari- 
son to the case of no interaction if a/kT ^ 1. This is op- 
posite to our findings. Quite probably, the disagreement 
stems from the under-relaxation of the initial (random) 
carrier configuration used in Ref. ; the relaxation pro- 
cess is pretty slow for interacting carriers if a/kT 3> 1. 
We cannot make more detailed comparison because typi- 
cal relaxation times are not provided in Ref. (in fact, 
even the strength of carrier repulsion is not provided). 
Our data indicates, for example, that for a/kT — 4 re- 
laxation is not completely over even for t/to = 1 x 10^ 
(see Fig. [3]) ; at that time carrier has already travelled in 
the field direction the distance of ~ 4 x lO'^a. 

Remarkable feature of Fig. [3] is a universality of the 
late relaxation stage. Very early relaxation is different 
for the initial random distribution and minimal energy 
distribution (where every carrier was placed at the site 
where the total energy, provided by the intrinsic disor- 
der and all previously added carriers, has a minimum), 
but after t/t^ ~ 10 relaxation curves merge into a single 
curve. 

Mobility field dependence 

It is very well known that the GDM and DG model 
provide very different dependences of the mobility on ap- 
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FIG. 3: Relaxation of the mean carrier velocity in the GDM 
for cr/kT — 4, eaE/a = 1, and p = 0.0625 (squares). Empty 
squares show relaxation for random initial locations, while 
filled squares show relaxation for the case, when initial loca- 
tions have been taken at the minimal energy positions. Circles 
shoe the realaxation for higher temperature a jkT = 3.33. 

plied electric field in the case of low carrier density. In 
the GDM log/i cx and in DG log/x (x @. These 
dependences remains valid for moderate -p too (see Fig. 

Significant difference between models could be found 
in the dependence of the slope of il{E) on p. General 
tendency for the DG model is that transformation of the 
mobility curve with the increase of p roughly resembles 
the corresponding transformation of the curve with the 
increase of T (see, for example. Fig. 1 in Ref. [ISj]): with 
the increase of p mobility becomes greater and the slope 
of the mobility field dependence becomes smaller. This is 
not the case for the GDM: here only mobility curve moves 
upward but the slope remains approximately constant. 

This difference could be easily understood. Field de- 
pendence of in the GDM is governed by the carrier es- 
cape from deep states to the nearest transport sites usu- 
ally having much higher energy and could be described 
by the shift AC/ = eER of positions of transport levels 
in applied field E 

log^oc-^, (6) 

where i? ~ a is a typical distance between transport sites. 
Random charge distribution provides a smooth random 
energy landscape superimposed on the intrinsic disorder, 
but typical additional variation of energy at the scale a 
is negligible for small p. Hence, estimation ^ remains 
valid and the slope of the mobility curve does not depend 
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FIG. 4: Mobility field dependence for Uc = 5, kTja = 0.3, 
and different values of p: 0.0016, 0.008, 0.016, and 0.032, from 
the bottom curve upward, correspondingly. Figure (a) shows 
simulation data for the DG model, and the figure (b) shows 
data for the GDM, correspondingly. Here /io ~ ea^To/a. If 
a = 1 nm and cr = 0.1 eV, then eaE ^ 1 ior E — 1 x 10® 
V/cm. 



on p. 

Situation in the DG model is quite different. Here mo- 
bility field dependence is governed by the carrier escape 
from wide valleys, and the size of critical valleys, capa- 
ble of keeping a carrier for the longest time, depends on 



5 



° -2 



Uc 


JU ^ . 

■20'^'^ 




5 g::::^ 
1 

0.3"^^ 













TOF 





0.5 



1 

(eaE/a) 



1.5 



1/2 



FIG. 5: Mobility field dependence for the DG model for dif- 
ferent values of Uc, indicated at the curves. Other parameters 
are: kT/a = 0.3 and p = 0.008. The very bottom curve cor- 
responds to the data of a conceivable TOF experiment with 
p^O. 
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Estimation ([6]) should be replaced by 
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and this result explains the origin of the famous Poole- 
Frenkel mobility field dependence observed in organic 
materials for the case of low density of carriers. If we in- 
crease the density of carriers, then at first they fill these 
critical traps, because the release time is maximal here. 
Hence, transport of more mobile carriers is governed by 
valleys (clusters) with the size that differs from i?cr- We 
can safely assume that the relevant traps has R < i?cr, 
because the distribution of clusters on size R decays with 
R. Hence, effective size R in Eq. ([5]) is smaller than i?cr, 
and, naturally, the slope of the mobility field dependence 
should be smaller than the slope of the corresponding 
curve in the case of p — > (and should decrease with the 
increase of p), exactly as it was observed in simulation. 

Fig. [5] illustrates transformation of the mobility curve 
with the increase of the strength of interaction. Again, 
effective interaction strength Uc affects mobility in the 
same way as temperature. It is very interesting to com- 
pare three curves in this figure: curves for Uc — 0, 
Uc — 5, and for a conceivable time-of-flight (TOF) ex- 
periment, i.e. for p — > 0. The first curve exactly corre- 
sponds to the results of Refs. P, Q, where only filling 



of deep states by non-interacting carriers has been taken 
into account, while the second one is much more close 
to the real situation. Evidently, when we consider the 
effect of carrier density on the mobility, then the effect of 
charge-charge interactions is no less important than the 
trivial filling of deep states. 



HOW TO COMPARE WITH EXPERIMENT? 



One can suggest that carrier transport in organic field- 
effect transistors (OFETs) should be a natural choice 
for comparison of the simulation results with experiment 
15|. Estimation of the carrier density in OFETs 
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show that the density as high as 3 x 10^^ cm ^ could be 
achieved [l^, that for a ~ 1 nm corresponds to p « 0.03. 
Experimental data for the particular OFET should be 
compared with the TOF data for a sandwich device hav- 
ing transport layer of the same material; in this way we 
could measure transport characteristics (e.g., cr), rele- 
vant for the intrinsic disorder in the material. Quite fre- 
quently, OFETs demonstrate mobilities much higher that 
the mobilities measured in TOF experiments, and usu- 
ally mobility increases with the increase of p [3] ■ This 
fact is in general agreement with the model studied in 
this paper. 

However, careful analysis reveals much more compli- 
cated situation. Indeed, in many aspects OFETs are 
very far away from the model, considered in the cur- 
rent study. First of all, in OFETs transport occurs in a 
thin layer, close to the gate insulator. Quite probably, 
especially in polymer devices, structure of this layer dif- 
fers from the structure of the same material in the bulk 
(polymer chains could be arranged in a special way at the 
gate insulator surface). This arrangement could provide 
more ordered structure with less degree of energetic dis- 
order, thus mobility should be enhanced near the inter- 
face, but accumulation of surface defects and impurities 
at the interface could lead to the decrease of the mobility. 
Next, there is a clear indication that the roughness of the 
organic semiconductor/dielectric interface affects carrier 
mobility At last, the very nature of a gate dielec- 

tric (specifically, its polarity) affects carrier mobility in 
OFETs, because a random orientation of polar groups in 
the vicinity of a transport layer induces an additional en- 
ergetic disorder in semiconductor 18,1^. In short, there 
are a lot of reasons to believe that transport properties of 
OFETs are too complicated to be directly compared with 
the results of this study. We can only state that a sig- 
nificant increase of the carrier mobility with the increase 
of carrier density in carefully manufactured OFETs does 
not contradict the results of our study. 
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CONCLUSION 

In conclusion, it was found that charge-charge interac- 
tion significantly affects the carrier drift mobility and this 
effect could not be described by the mobility dependence 
on average electric field in the model of non-interacting 
carriers. Effect of the charge-charge interaction is no less 
important for the true description of the carrier mobility 
dependence on carrier density than the previously con- 
sidered effect of simple filling of deep states. 
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